Formulations of the Einstein equations for numerical simulations * 
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We review recent efforts to re-formulate the Einstein equations for fully relativistic numerical 
simulations. The so-called numerical relativity is a promising research field matching with ongo- 
ing gravitational wave observations. In order to complete a long-term and accurate simulations of 
binary compact objects, people seek a robust set of equations against the violation of constraints. 
Many trials have revealed that mathematically equivalent sets of evolution equations show different 
numerical stability in free evolution schemes. In this article, we overview the efforts of the com- 
munity, categorizing them into three directions: (1) modifications of the standard Arnowitt-Deser- 
Misner equations initiated by the Kyoto group (the so-called Baumgarte-Shapiro-Shibata-Nakamura 
equations), (2) rewriting the evolution equations in a hyperbolic form, and (3) construction of an 
"asymptotically constrained" system. We then introduce our series of works that tries to explain 
these evolution behaviors in a unified way using eigenvalue analysis of the constraint propagation 
equations. The modifications of (or adjustments to) the evolution equations change the character of 
constraint propagation, and several particular adjustments using constraints are expected to damp 
the constraint- violating modes. We show several set of adjusted ADM/BSSN equations, together 
with their numerical demonstrations. 

PACS numbers: 04.20.-q, 04.20.Cv, 04.25.D- 
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I. INTRODUCTION 



Overview 



The theory of general relativity describes the nature of 
the strong gravitational field. The Einstein equation pre- 
dicts quite unexpected phenomena such as gravitational 
collapse, gravitational waves, the expanding universe and 
so on, which are all attractive not only for researchers but 
also for the public. The Einstein equation consists of 10 
partial differential equations (elliptic and hyperbolic) for 
10 metric components, and it is not easy to solve them for 
any particular situation. Over the decades, people have 
tried to study the general-relativistic world by finding 
its exact solutions, by developing approximation meth- 
ods, or by simplifying the situations. While "The Exact 
Solution" book [70| says there were more than 4000 pub- 
lications on exact solutions between 1980 and 2000, di- 
rect numerical integration of the Einstein equations can 
be said to be the most robust way to study the strong 
gravitational field. This research field is often called "nu- 
merical relativity" . 

With the purpose of the predictions of precise grav- 
itational waveforms from coalescences of the binary 
neutron-stars and/or black-holes, numerical relativity 
has been developed for the past three decades. The diffi- 
culty of numerical integrations of the Einstein equations 
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arises from its mathematical complexity of the equations, 
physical difficulty of singularity treatments, and from 
high-level requirements for computational skills and tech- 
nology. 

In 2005-2006, several groups independently announced 
that the simulations of the inspiral black-hole binary 
merger are available [3, [1^, [3l|, [H, [6^ . There are many 
implements for their successes, such as gauge conditions, 
coordinate selections, boundary treatments, singularity 
treatments, numerical discretization, and mesh refine- 
ments, together with the re-formulation of the Einstein 
equations which we will discuss here. More general and 
recent introductions to numerical relativity are available, 
e^. by Baumgarte-Shapiro (2003) [ll|, Alcubierre (2004) 
Pretorius (2007) 55], and Bruegmann (2008) 0]. 

The purpose of this article is to review the formula- 
tion problem in numerical relativity. This is one of the 
essential issues to achieve the long-term stable and accu- 
rate simulations of binary compact objects. Mathemati- 
cally equivalent sets of evolution equations show different 
numerical stability in free evolution schemes. This had 
been the mystery for long time between the relativists, 
and many proposals and trials were reported. After we 
review the problem from such a historical viewpoint, we 
will explain our systematic understanding using the con- 
straint propagation equations; the evolution equations 
of the constraints which is supposed to be satisfied all 
through the time evolutions. 

The most numerical relativity groups today uses the 
so-called BSSN (Baumgarte-Shapiro-Shibata-Nakamura) 
equations, that is one of the modified form of the ADM 
(Arnowitt-Deser-Misner) equations. We try to explain 
why these differences appear and also predict that more 
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robust sets of equations exist together with actual nu- 
merical demonstrations. 



B. Formulation Problem in Numerical Relativity 

There are several different approaches to simulating 
the Einstein equations. Among them the most robust 
way, that we target in this article, is to apply 3+1 (space 
+ time) decomposition of space-time. This formulation 
was first given by Arnowitt, Deser and Misner (ADM) 
[T^ l (we call the original ADM system, hereafter) with 
the purpose of constructing a canonical formulation of 
the Einstein equations to seek the quantum nature of 
space-time. In late 70s, when the numerical relativity 
started, this ADM formulation was introduced by Smarr 
and York [gl, [t^ in a slightly different notations which 
is taken as the standard formulation between numerical 
relativists (so that we call the standard ADM system, 
hereafter) . 

The ADM formulation divides the Einstein equations 
into the constraint equations and the evolution equations 
apparently, like the Maxwell equations. Since the set of 
ADM equations form a first-class system, if we solved 
two constraint equations, the Hamiltonian (or energy) 
constraint and the momentum constraint equations for 
the initial data, then the evolution equations theoreti- 
cally guarantees the evolved data satisfy the constraint 
equations. This free- evolution approach is also the stan- 
dard in numerical relativity. This is because solving the 
constraints (non-linear elliptic equations) is numerically 
expensive, and because the free-evolution allows us to 
monitor the accuracy of numerical evolution using the 
constraint equations. 

Up to the middle 90s, the ADM numerical relativity 
appealed great successes. For example, the formation 
of naked singularity from collisionless particles [sOl shows 
the unknown behavior of the strong gravity; the discovery 
of the critical behavior for a black-hole formation [2g| 
open-doors the understanding of phase-transition nature 
in general relativity; the black-hole horizon dynamics 
realized the theoretical predictions. 

Nevertheless, when people try to make the long-term 
simulations such as coalescences of neutron-star binaries 
and/or black- hole binaries for calculating gravitational- 
wave form, numerical simulations were often interrupted 
by unexplained blow-ups (Figure [1]). This was thought 
due to the lack of resolution, or inappropriate gauge 
choice, or the particular numerical scheme which was ap- 
plied. However, with the accumulation of experiences, 
people have noticed the importance of the formulation of 
the evolution equations, since there are apparent differ- 
ences in numerical stability although the equations are 
mathematically equivalent. 

At this moment, there are three major ways to obtain 
longer time evolutions, which we describe in the next 
section. Of course, the ideas, procedures, and problems 
are mingled with each other. The purpose of this article 



Blowup 




FIG. 1; Origin of the problem for numerical relativists: Nu- 
merical evolutions depart from the constraint surface. 



is to review all three approaches and to introduce our 
idea to view them in a unified way. The author wrote a 
detail review of this topic in 2002 [HH], and the present 
article includes an update in brief style. 

The word stability is used quite different ways in the 
community. 

• We mean by numerical stability a numerical simu- 
lation which continues without any blow-ups and 
in which data remains on the constrained surface. 

• Mathematical stability is defined in terms of the 
well-posedness in the theory of partial differential 
equations, such that the norm of the variables is 
bounded by the initial data. See eq. (^5)1 and 
around. 

• For numerical treatments, there is also another no- 
tion of stability, the stability of finite differencing 
schemes. This means that numerical errors (trun- 
cation, round-off, etc) are not growing by evolution. 
The evaluation is obtained using von Neumann's 
analysis. Lax's equivalence theorem says that if a 
numerical scheme is consistent (converging to the 
original equations in its continuum limit) and sta- 
ble (no error growing), then the simulation repre- 
sents the right (converging) solution. See [131 for 
the Einstein equations. 

We follow the notations of that of MTWjlOl. We use 
fi, v ~ 0, ■ ■ ■ ,3 and i, j = 1, • • • ,3 as space-time indices. 
The unit c = 1 is applied. The discussion is mostly to 
the vacuum space-time, but the inclusion of matter is 
straightforward. 



3 



surface normal line 



shift vectoiN' 




2. The standard ADM formulation 



t = constant hyperstirf ace 



FIG. 2: Concept of time evolution of space-time: foliations 
of 3-dimensional hypersurface. The lapse and shift functions 
are often denoted a or A'', and or respectively. 



II. THE STANDARD WAY AND THE THREE 
OTHER ROADS 

A. Strategy 0: The ADM formulation 

1. The original ADM formulation 

The Arnowitt-Deser-Misner (ADM) formulation 
gave the fundamental idea of time evolution of space and 
time: such as foliations of 3-dimcnsional hypersurface 
(Figured]). 

The story begins by decomposing 4-dimensional space- 
time into 3 plus 1. The metric is expressed by 



^dt^ + 7y {dx' + I3'dt){dx^ -t- 13^ dt), (1) 



where a and (3j are defined as a 



1 



and 



f3j = goj, and called the lapse function and shift vec- 
tor, respectively. The projection operator or the intrin- 
sic 3- metric gij is defined as 7^,^ = g^,, + n^n^, where 
= (-a, 0,0,0) [and = gf^^n^ = {l/a,-(3'/a)] is 
the unit normal vector of the spacelike hypersurface, E 
(see Figure [2]). By introducing the extrinsic curvature. 



(2) 



and using the Gauss-Codacci relation, the Hamiltonian 
density of the Einstein equations can be written as 

Hgr = 7^''i^J - (3) 

where 

£ = = a^l^^^R -K^ + K,jK''] , (4) 

where tt*-' is the canonically conjugate momentum to 7^ , 

-^{K^^ - (5) 



dC 

dii3 



omitting the boundary terms. The variation of H-gr with 
respect to a and f3i yields the constraints, and the dynam- 
ical equations are given by 7^^ = ^j^^ and -k^^ — — ^^^^ 



Kij was 



In the version of Smarr and York[68|, |7l 
used as a fundamental variable instead of the conjugate 
momentum tt*-' . The set of equation is summarized as 
follows: 

□ The Standard ADM formulation fs^ . It^ / 
The fundamental dynamical variables are {'jij,Kij), 
the three-metric and extrinsic curvature. The three- 
hypersurface E is foliated with gauge functions, (a,/3*), 
the lapse and shift vector. 



The evolution equations: 



(6) 



-r oiKKij — 2aKikK j — DiDja 

where K — K^i, and '■^-'^y and denote 
three-dimensional Ricci curvature, and a covariant 
derivative on the three-surface, respectively. 



• Constraint equations: 



j^ADM 
^ADM 



+ K'^ - KuK'^ w 0, 



D,K » 0, 



(8) 
(9) 



where ^^'^R =(3) Jii^. these are called the Hamilto- 
nian (or energy) and momentum constraint equa- 
tions, respectively. □ 

The formulation has 12 first-order dynamical variables 
(7ij , Kij)^ with 4 freedom of gauge choice (a, [3i) and with 
4 constraint equations, ([8]) and The rest freedom 

expresses 2 modes of gravitational waves. 

We remark that there is one replacement in ([7|) using 
([8|) in the process of conversion from the original ADM 
to the standard ADM equations. This is the key issue in 
the later discussion, and we shall be back this point in 
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The constraint propagation equations, which are the 
time evolution equations of the Hamiltonian constraint 
([8|) and the momentum constraints ([9]), can be written 
as follows: 

OConstraint Propagations of the Standard ADM: 

dtH = iS^djH) + 2aKn - 2aY^ (d^Mj) 

-47^^(9,a)X„ (10) 
^tM^ = -{^/2)ai^^n)-{^,a)n + (3Ji^JM^) + aKM^ 

-P''i''{d^jik)Mj + {dAh'''Mj. a (11) 

From these equations, we know that if the constraints 
are satisfied on the initial slice E, then the constraints 
are satisfied throughout evolution. The normal numerical 
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FIG. 3: Chronological table of formulations and their numerical tests (~ 2001). Boxed ones are of proposals of formulation, 
circled ones are related numerical experiments. Please refer Table 1 in [651 ] for each references. 



scheme is to solve the elliptic constraints for preparing 
the initial data, and to apply the free evolution (solving 
only the evolution equations). The constraints are used 
to monitor the accuracy of simulations. 

The ADM formulation was the standard formulation 
for numerical relativity up to the middle 90s. Numerous 
successful simulations were obtained for the problems of 
gravitational collapse, critical behavior, cosmology, and 
so on. However, stability problems have arisen for the 
simulations such as the gravitational radiation from com- 
pact binary coalescence, because the models require quite 
a long-term time evolution. 

The origin of the problem was that the above state- 
ment in Italics is true in principle, but is not always true 
in numerical applications. A long history of trial and er- 
ror began in the early 90s. From the next subsection we 
shall look back of them by summarizing "three strate- 
gies" . We then unify these three roads as "adjusted sys- 
tems" , and as its by-product we show that the standard 
ADM equations has a constraint violating mode in its 
constraint propagation equations even for a single black- 
hole (Schwarzschild) spacetime [64]. Figure [3] and H] are 
chronological maps of the researches. 



B. Strategy 1: Modified ADM formulation by 
Nakamura et al (The BSSN formulation) 

Up to now, the most widely used formulation for 
large scale numerical simulations is a modified ADM 
system, which is now often cited as the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) formulation. This 



re-formulation was first introduced by Nakamura et al. 
[53 , m, [6l[. The usefulness of this re- formulation was 
re- introduced by Baumgarte and Shapiro |15|, then was 
confirmed by other groups to show a long-term stable 
numerical evolution [j, Q . 



1. Basic variables and equations 

The widely used notation [l^ introduces the variables 
{(p,jij,K,Aij,r'') instead of {'yij,Kij), where 
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log(det7. 



Uj = e-^^{K,, - ( 1/3)7., X), 



(12) 

(13) 
(14) 
(15) 
(16) 



The new variable is introduced in order to calcu- 
late Ricci curvature more accurately. In BSSN formu- 
lation, Ricci curvature is not calculated as Rfj^^^ = 

^k^ij~^i^kj~^^ij^fk~^kj^H^ but aS Rfj^^^ = Rfj+Rij, 

where the first term includes the conformal factor ip while 
the second term does not. These are approximately 
equivalent, but Rfj^^'^ does have wave operator appar- 
ently in the flat background limit, so that we can expect 
more natural wave propagation behavior. 

Additionally, the BSSN requires us to impose the con- 
formal factor as 7(:= det^ij) = 1, during evolution. This 



5 



2001 



2005 




CUTB^Rochester) ^ ^ 

7^==:^^ (Parma^ 



BSSN is "well-posed" ? 
CSarbach / Gundlach . . . ) 



asymptotically constrained / 
constraint damping 



Qiana et.al.) 

'''>k JZ4-lQmbdQl ^s 
_j_ CGundlach-Calabrese) 

»— 'C ^nko'i -Yone do^y (^^^Pretoriur^ 

/ 




CYoneda-Shinkai) 



CShinkai-Yoneda)»~«CShinkai 



FIG. 4: Chronological table of formulations and their numerical tests (2001 
circled ones are related numerical experiments. 



Boxed ones are of proposals of formulation, 



is a kind of definition, but can also be treated as a con- 
straint. 

□ The BSSN formulation fM, \M, \M, \M]: 

The fundamental dynamical variables are 

The three-hypersurface E is foliated with gauge func- 
tions, (a,/3'), the lapse and shift vector. 

• The evolution equations: 

= -(l/6)aX+(l/6)/3X5.^) + (9./3'), (17) 

~{2/3)j,,idkP'') + P'{dk%), (18) 

- -e-4^(Ai?j«)^^ + e-4^a(i?,f ^^)^^ 

+aKA,j - 2ai,fci^- + (9,/3'=)ifc, -f (9j/3'=)ife, 
-(2/3)(afe/3'=)iy + P^dkA,,), (20) 



-2(9,a)i»^- + 2a(f;.,i'=-'- - (2/3)r-'(a,X) 
+6i^^"(9,(^)) - d,{(3''{dkf') - 7"^(afc/30 
-7'=Xafe/3^) + (2/3)^^(9^/3'=)). (21) 



• Constraint equations: 

^BSSN ^ ^BSSN ^ _ j^^^j^ij ^ 
j^BSSN _ j^ADM 

= r-7^''=f;.;,, 

A = A,jf\ 

5 = 7-1. □ 



(22) 
(23) 
(24) 

(25) 
(26) 



and (|23p are the Hamiltonian and momentum con- 
straints (the "kinematic" constraints), while the latter 
three are "algebraic" constraints due to the requirements 
of BSSN variables. 



2. Remarks, Pros and Cons 

Why BSSN is better than the standard ADM? To- 
gether with numerical comparisons with the standard 
ADM easeful, this question has been studied by many 
groups using different approaches. 

• Using numerical test evolutions, Alcubierre et al. 
pi] found that the essential improvement is in the 
process of replacing terms by the momentum con- 
straints. They also pointed out that the eigenvalues 
of BSSN evolution equations have fewer "zero eigen- 
values" than those of ADM, and they conjectured 
that the instability might be caused by these "zero 
eigenvalues" . 

• Miller [4^ reported that BSSN has a wider range of 
parameters that gives us stable evolutions in von 
Neumann's stability analysis. 

• An effort was made to understand the advantage 
of BSSN from the point of hyperbolization of the 
equations in its linearized limit 0, [s^l, or with a 
particular combination of slicing conditions plus 
auxiliary variables [isj. If we define the 2nd order 
symmetric hyperbolic form, then the principal part 
of BSSN can be one of them[41]. 
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As we discussed in [tJ, the stability of the BSSN for- 
mulation is due not only to the introductions of new vari- 
ables, but also to the replacement of terms in the evo- 
lution equations using the constraints. Further, we can 
show several additional adjustments to the BSSN equa- 
tions which give us more stable numerical simulations. 
We will devote jjllll for the fundamental idea. 

The current binary black-hole simulations apply the 
BSSN formulations with several implementations. For 
example, 

tip-1 Alcubierre et al. reported that trace-out Aij 
technique at every time-step helps the stability. 

tip-2 Campanelli et al. reported that in their code f ' 
is replaced by —djj^^ where it is not differentiated. 

tip-3 Baker et al. [l^ modified F'-equation ([2T|l as sug- 
gested by Yo et al. [73| . 

These technical tips are again explained using the con- 
straint propagation analysis as we will come back in 

micii 

These studies provide sort of evidences regarding the 
advantage of BSSN, while it is also shown an example 
of an ill-posed solution in BSSN (as well in ADM) by 
Frittelli and Gomez [s^]. Recently, the popular combi- 
nation, BSSN with Bona-Masso type slicing condition, is 
investigated in particular. Among then, Garfinkle et al. 
[39I I speculated that the reason of gauge shocks are miss- 
ing in the current 3-dimensional black-hole simulations is 
simply because the lack of resolution. 



C. Strategy 2: Hyperbolic re-formulations 

1. Definitions, properties, mattiematical bacl:grounds 

The second effort to re-formulate the Einstein equa- 
tions is to make the evolution equations reveal a first- 
order hyperbolic form explicitly. This is motivated by the 
expectation that the symmetric hyperbolic system has 
well-posed properties in its Cauchy treatment in many 
systems and also that the boundary treatment can be 
improved if we know the characteristic speed of the sys- 
tem. 

OHyperbolic formulations 

We say that the system is a first-order (quasi-linear) 
partial differential equation system, if a certain set of 
(complex- valued) variables Ua (a = 1, • • • ,n) forms 

dtUo,^M'^c.{u)diUi3+J^c.{u), (27) 

where M (the characteristic matrix) and M are functions 
of u but do not include any derivatives of u. Further we 
say the system is 

• a weakly hyperbolic system, if all the eigenvalues of 
the characteristic matrix are real. 



• a strongly hyperbolic system (or a diagonalizable / 
symmetrizable hyperbolic system) , if the character- 
istic matrix is diagonalizable (has a complete set of 
eigenvectors) and has all real eigenvalues. 

• a symmetric hyperbolic system, if the characteristic 
matrix is a Hermitian matrix. □ 

Writing the system in a hyperbolic form is a quite use- 
ful step in proving that the system is well-posed. The 
mathematical well-posedness of the system means (1°) 
local existence (of at least one solution u), (2°) unique- 
ness (i.e., at most solutions), and (3°) stability (or con- 
tinuous dependence of solutions {u} on the Cauchy data) 
of the solutions. The resultant statement expresses the 
existence of the energy inequality on its norm, 

IKi)||<e"nKi = 0)||, 
where < t < i, a = const. (28) 

This indicates that the norm of u{t) is bounded by a 
certain function and the initial norm. Remark that this 
mathematical bounds does not mean that the norm u{t) 
decreases along the time evolution. 

The inclusion relation of the hyperbolicities is, 

symmetric hyperbolic C strongly hyperbolic 

C weakly hyperbohc. (29) 

The Cauchy problem under weak hyperbolicity is not, in 
general, C°° well-posed. At the strongly hyperbolic level, 
we can prove the finiteness of the energy norm if the char- 
acteristic matrix is independent of u (cf [Tlj). that is one 
step definitely advanced over a weakly hyperbolic form. 
Similarly, the well-posedness of the symmetric hyperbolic 
is guaranteed if the characteristic matrix is independent 
of u, while if it depends on u we have only limited proofs 
for the well-posedness. 

From the point of numerical applications, to hyper- 
bolize the evolution equations is quite attractive, not only 
for its mathematically well-posed features. The expected 
additional advantages are the following. 

(a) It is well known that a certain flux conservative 
hyperbolic system is taken as an essential formula- 
tion in the computational Newtonian hydrodynam- 
ics when we control shock wave formations due to 
matter. 

(b) The characteristic speed (eigenvalues of the prin- 
cipal matrix) is supposed to be the propagation 
speed of the information in that system. There- 
fore it is naturally imagined that these magnitudes 
are equivalent to the physical information speed of 
the model to be simulated. 

(c) The existence of the characteristic speed of the sys- 
tem is expected to give us an improved treatment 
of the numerical boundary, and/or to give us a 
new well-defined Cauchy problem within a finite re- 
gion (the so-called initial boundary value problem; 
IB VP). 
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These statements sound reasonable, but have not yet 
been generally confirmed in actual numerical simulations 
in general relativity. 



2. Hyperbolic formulations of the Einstein equations 

Most physical systems can be expressed as symmetric 
hyperbolic systems. In order to prove that the Einstein's 
theory is a well-posed system, to hyperbolize the Einstein 
equations is a long-standing research area in mathemat- 
ical relativity. 

The standard ADM system does not form a first order 
hyperbolic system. This can be seen immediately from 
the fact that the ADM evolution equation ([7]) has Ricci 
curvature in RHS. This is also the common fact to the 
BSSN formulation. 

So far, several first order hyperbolic systems of the Ein- 
stein equation have been proposed. In constructing hy- 
perbolic systems, the essential procedures are (1°) to in- 
troduce new variables, normally the spatially derivatived 
metric, (2°) to adjust equations using constraints. Occa- 
sionally, (3°) to restrict the gauge conditions, and/or (4°) 
to rescale some variables. Due to process (1°), the num- 
ber of fundamental dynamical variables is always larger 
than that of ADM. 

Due to the limitation of space, we can only list several 
hyperbolic systems of the Einstein equations: 

• The Bona-Masso formulation 

• The Einstein-Ricci system [l|, / Einstein- 
Bianchi system Q 

• The Einstein-Christoffel system [l3| 

• The Ashtekar formulation [l^, (tI] 

• The Frittelli-Reula formulation [H, [Zll 

• The Conformal Field equations [s^ 

• The Bardeen-Buchman system [l3| 

• The Kidder-Schccl-Teukolsky (KST) formulation 
M 

• The Alekseenko- Arnold system Q 

• The general-covariant Z4 system [l^ 

• The Nagy-Ortiz-Reula (NOR) formulation [5l| 

• The Weyl system [H, IH 

Note that there is no apparent differences between the 
word 'formulation' and 'system' here. 



3. Numerical tests 

When we discuss hyperbolic systems in the context 
of numerical stability, the following questions should be 
considered: 

Q From the point of the set of evolution equations, 
does hyperbolization actually contribute to numer- 
ical accuracy and stability? Under what condi- 
tions/situations will the advantages of hyperbolic 
formulation be observed? 

Unfortunately, we do not have conclusive answers to 
these questions, but many experiences are being accumu- 
lated. Several earlier numerical comparisons reported the 
stability of hyperbolic formulations [13, [2l|, Ha, [HI . But 
we have to remember that this statement went against 
the standard ADM formulation. 

These partial numerical successes encouraged the com- 
munity to formulate various hyperbolic systems. How- 
ever, several numerical experiments also indicate that 
this direction is not a complete success. 

• Above earlier numerical successes were also termi- 
nated with blow-ups. 

• If the gauge functions are evolved according to the 
hyperbolic equations, then their finite propagation 
speeds may cause pathological shock formations in 
simulations [1, [1|. 

• There are no drastic differences in the evolu- 
tion properties between hyperbolic systems (weakly, 
strongly and symmetric hyperbolicity) by system- 
atic numerical studies by Hern [42l | based on 
Frittelli-Reula formulation [33| , and b y th e authors 
[63j based on Ashtekar's formulation fisl. [t^. 

• Proposed symmetric hyperbolic systems were not 
always the best ones for numerical evolution. Peo- 
ple are normally still required to re-formulate them 
for suitable evolution. Such efforts are seen in the 
applications of the Einstein-Ricci system the 
Einstein-Christoffel system flTj, and so on. 

• If we can erase the non-principal part by suit- 
able re-definition of variables (as is in the KST 
formulation) [4^ , then we can see the growth of 
energy norm (j28p in numerical evolution as theo- 
retically predicted [H, \^ . We then see a certain 
differences in the long-term convergence features 
between weakly and strongly hyperbolic systems. 

Of course, these statements only casted on a particular 
formulation, and therefore we have to be careful not to 
over-emphasize the results. 
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4- Remarks 

In order to figure out the reasons for the above objec- 
tions, it is worth stating the following cautions: 

(a) Rigorous mathematical proofs of well-posedness of 
PDE are mostly for simple symmetric or strongly 
hyperbolic systems. If the matrix components or 
coefficients depend on dynamical variables (as in all 
any versions of hyperbolized Einstein equations), 
almost nothing was proved in more general situa- 
tions. 

(b) The statement of "stability" in the discussion of 
well-posedness refers to the bounded growth of the 
norm, ([28]) . and does not indicate a decay of the 
norm in time evolution. 

(c) The discussion of hyperbolicity only uses the char- 
acteristic part of the evolution equations, and ig- 
nores the rest. 



Blow up 




FIG. 5: The idea of the "Asymptotically Constrained Sys- 
tem" . 



We think the origin of confusion in the community re- 
sults from over-expectation on the above issues. Mostly, 
point (c) is the biggest problem. The above numerical 



claims from Ashtekar(63|, |75l and Frittelli-Reula [43 for- 
mulations were mostly due to the contribution (or inter- 
position) of non-principal parts in evolution. Regarding 
this issue, the KST formulation finally opens the door. 
KST's "kinematic" parameters enable us to reduce the 
non-principal part, so that numerical experiments are 
hopefully expected to represent predicted evolution fea- 
tures from PDE theories. At this moment, the agreement 
between numerical behavior and theoretical prediction is 
not yet perfect but close [isj . 

If further studies reveal the direct correspondences be- 
tween theories and numerical results, then the direction 
of hyperbolization will remain as the essential approach 
in numerical relativity, and the related IBVP researches 
[H, il m HE [zH will become a main research subject 
in the future. Meanwhile, it will be useful if we have an 
alternative procedure to predict stability including the 
effects of the non-principal parts of the equations. Our 
proposal of adjusted system in the next subsection may 
be one of them. 



1, The "X-system" 

Brodbeck et al. [l^l proposed a system which has ad- 
ditional variables A that obey artificial dissipative equa- 
tions. The variable As are supposed to indicate the vi- 
olation of constraints and the target of the system is to 
get A = as its attractor. 
DThe "X-system" (Brodbeck et al.) fMj : 
For a symmetric hyperbolic system, add additional vari- 
ables A and artificial force to reduce the violation of con- 
straints. The procedure is as follows: 

1. Prepare a symmetric hyperbolic evolution system 



dtu = MdiU + N 



(30) 



2. Introduce A as an indicator of violation of con- 
straint which obeys dissipative eqs. of motion 

dtX = aC-(3X,{a^0,(3>0) (31) 

3. Take a set of (u, A) as dynamical variables 

(32) 



dt 



A 0\„ (u 



D. Strategy 3: Asymptotically constrained systems 

The third strategy is to construct a robust system 
against the violation of constraints, such that the con- 
straint surface is an attractor (Figure [5]) . The idea was 
first proposed as "A-system" by Brodbeck et al. , and 
then developed in more general situations as "adjusted 
system" by the authors [75^. 



4. Modify evolution eqs so as to form a symmetric 
hyperbolic system 

*(°)-(fo)<a) ° ('^) 

Since the total system is designed to have symmetric hy- 
perbolicity, the evolution is supposed to be unique. Brod- 
beck et al. showed analytically that such a decay of As 
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can be seen for sufHciently small A(> 0) with a choice of 
appropriate combinations of as and /?s. 

Brodbeck et al. presented a set of equations based on 
Frittelli-Reula's symmetric hyperbolic formulation [ss! ]. 
The version of Ashtekar's variables was presented by the 
authors [12] for controlling the constraints or reality con- 
ditions or both. The numerical tests of both the Maxwell- 
A-system and the Ashtekar-A-system were performed [TSj , 
and confirmed to work as expected. The A-system version 
of the general-covariant Z4 system [l^ is also presented 
(40| . Pretorius Q applied this "constraint-damping" 
idea in his harmonic system to perform his binary black- 
hole merger simulations. 

Although it is questionable whether the recovered so- 
lution is true evolution or not j67| , we think the idea is 
quite attractive. To enforce the decay of errors in its 
initial perturbative stage seems the key to the next im- 
provements, which are also developed in the next section 
on "adjusted systems" . 

However, there is a high price to pay for constructing 
a A-system. The A-system can not be introduced gen- 
erally, because (i) the construction of A-system requires 
the original evolution equations to have a symmetric hy- 
perbolic form, which is quite restrictive for the Einstein 
equations, (ii) the final system requires many additional 
variables and we also need to evaluate all the constraint 
equations at every time step, which is a hard task in com- 
putation. Moreover, (iii) it is not clear that the A-system 
is robust enough for non-linear violation of constraints, 
or that A-system can control constraints which do not 
have any spatial differential terms. 

2. The "adjusted system" 

Next, we propose an alternative system which also tries 
to control the violation of constraint equations actively, 
which we named "adjusted system" . We think that this 
system is more practical and robust than the previous 
A-system. The essential is summarized as follows: 
□ The Adjusted system (procedures): 

1. Prepare a set of evolution eqs. 

dtu = Jd,u + K (34) 

2. Add constraints in RHS 

dtu = Jd,u + K+K^ (35) 

3. Choose the coefficient (or Lagrange multiplier) k 
so as to make the eigenvalues of the homogenized 
adjusted dtC equations negative real value or pure 
imaginary. 

dtC ^ Dd^C + EC (36) 
dtC = DdjC + EC +FdjC + GC D (37) 



The process of adjusting equations is a common tech- 
nique in other re- formulating efforts as we reviewed. 
However, we try to employ the evaluation process of con- 
straint amplification factors as an alternative guideline 
to hyperbolization of the system. We will explain these 
issues in the next section. 



III. A UNIFIED TREATMENT: ADJUSTED 
SYSTEM 

This section is devoted to present our idea of "asymp- 
totically constrained system" . Original references can be 
found in [6l[7i,[7i,[73. 

A. Procedures : Constraint propagation equations 
and Proposals 

Suppose we have a set of dynamical variables (a;* , t) , 
and their evolution equations, 

9t^i'^ = /(^i^a,^i^•••), (38) 

and the (first class) constraints, 

C"(w°,5iu",---)«0. (39) 

Note that we do not require ([38|) forms a first order hy- 
perbolic form. We propose to investigate the evolution 
equation of C" (constraint propagation), 

9tC"=g(C",9,C",---), (40) 

for predicting the violation behavior of constraints in 
time evolution. We do not mean to integrate ([3D]) nu- 
merically together with the original evolution equations 
psp . but mean to evaluate them analytically in advance 
in order to re- formulate the equations ([55]) . 

There may be two major analyses of (HO]); (a) the hy- 
perbolicity of (|40)) when (|40l) is a first order system, and 
(b) the eigenvalue analysis of the whole RHS in after 
a suitable homogenization. As we mentioned in ^H C 4| 
one of the problems in the hyperbolic analysis is that it 
only discusses the principal part of the system. Thus, we 
prefer to proceed the road (b). 
OConstraint Amplification Factors (CAFs): 
We propose to homogenize (|40l) by a Fourier transforma- 
tion, e.g. 

where C(x, t)'' = J C{k,t)P exp{ik ■ x)d^k, (41) 

then to analyze the set of eigenvalues, say As, of the co- 
efficient matrix, M^p, in ([1T|) . We call As the constraint 
amphfication factors (CAFs) of (gO]) . □ 

The CAFs predict the evolution of constraint viola- 
tions. We therefore can discuss the "distance" to the 
constraint surface using the "norm" or "compactness" of 
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the constraint violations (although we do not have exact 
definitions of these "• • • " words). 

The next conjecture seems to be quite useful to predict 
the evolution feature of constraints: 
OConjecture on CAFs: 

(A) If CAF has a negative real-part (the constraints are 
forced to be diminished), then we see more stable 
evolution than a system which has positive CAF. 

(B) If CAF has a non-zero imaginary-part (the con- 
straints are propagating away), then we sec more 
stable evolution than a system which has zero CAF. 
□ 

We found that the system becomes more stable when 
more As satisfy the above criteria. (The first observations 
were in the Maxwell and Ashtekar formulations [63.. .75,]). 

Actually, supporting mathematical proofs are available 
when we classify the fate of constraint propagations as 
follows. 

n\ Classification of Constraint Propagations: 
If we assume that avoiding the divergence of constraint 
norm is related to the numerical stability, the next clas- 
sification would be useful: 

(CI) Asymptotically constrained : All the constraints 
decay and converge to zero. 

This case can be obtained if and only if all the real 
parts of CAFs are negative. 

(C2) Asymptotically bounded : All the constraints are 
bounded at a certain value, (this includes the above 
asymptotically constrained case.) 
This case can be obtained if and only if (a) all the 
real parts of CAFs are not positive and the con- 
straint propagation matrix M" p is diagonalizable, 
or (b) all the real parts of CAFs are not positive 
and the real part of the degenerated CAFs is not 
zero. 

(C3) Diverge : At least one constraint will diverge. 
□ 

The details are shown in [t^. A practical procedure for 
this classification is drawn in Figure [S] 

The above features of the constraint propagation, 
will differ when we modify the original evolution equa- 
tions. Suppose we add (adjust) the evolution equations 
using constraints 

dtu^' = /(u^ a,u^ • • • ) + F(c", • • • ), (42) 

then (pOj) will also be modified as 

dtC"' = .9(C",5iC",---) + G(C",9^C",---)- (43) 

Therefore, the problem is how to adjust the evolution 
equations so that their constraint propagations satisfy 
the above criteria as much as possible. 



Ql: Is there a CAP which real part is positive? 

NO / YES ». Diverge 



Q2: Are all the real parts of CAFs negative? 

NO / YES 



Asymptotically 
Constrained 



Q3: Is the constraint propagation matrix diagonalizable? 

NO / YES 



Asymptotically 
Bounded 



Q4: Is a real part of the degenerated CAFs is zero? 



YES / NO 



Asymptotically 
Bounded 



Q5: Is the associated Jordan matrix diagonal? 

NO / YES 

I 

Diverge 



Asymptotically 
Bounded 



FIG. 6: A flowchart to classify the constraint propagations. 



B. Applications 1: Adjusted ADM formulations 

1. Adjusted ADM equations 

Generally, we can write the adjustment terms to ^ 
and ([7]) using ([8]) and ([9]) by the following combinations 
(using up to the first derivatives of constraints for sim- 
plicity): 

OThe adjusted ADM formulation '61]: 

Modify the evolution eqs of {'-fij,Kij) by constraints H 

and ^Ai, i.e., 



dtKij 



+p\j{ykn) -f q''\jiVkMi), (44) 

III +R^Jn + S''^JMk 



,{VkH) + s'',^iVkM 



(45) 



where P, Q, R, S and p, g, r, s are multipliers. According 
to this adjustment, the constraint propagation equations 
are also modified as 

dtH (dni) + additional terms (46) 
dtMi = (HI]) -I- additional terms □ (47) 

We show two examples of adjustments here. Several 
others are shown in Table 3 of (64l |. 

1. The standard ADM vs. original ADM 

The first comparison is to show the differences 
between the standard ADM [t^ and the original 
ADM system fl^] (see ^11 A|) . In the notation of 
pil) and (gSl), the adjustment. 



(48) 



Rij = Kpa'jij, 



(and set the other multipliers zero) will distinguish 
two, where kf is a constant. Here /tj? = corre- 
sponds to the standard ADM (no adjustment), and 
Kp = —1/ A to the original ADM (without any ad- 
justment to the canonical formulation by ADM). 
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ADM (standard) 

Delweilei k=0,05. off <S 10"(-3) 
DelwBiler k=0 05, off @ 10"(-4) 
Detweiler k=0.05, off @ l0"(-5) 
Detweiler k=0.05 




time 



100 200 300 400 500 600 700 800 

time 



FIG. 7: Demonstration of numerical evolutions between adjusted ADM systems; especially the standard ADM system and 
Detweiler's modified ADM system. The L2 norm of the constraints Ji^^^' ig plotted in the function of time. The model 
is propagation of Teukolsky wave in the periodical 3-dimensional box. k is the parameter in Detweiler's adjustment [fci in 
ea. (l49p - (l52l) ]. with fixed-fc cases (left panel) and with fixed-and-turnoff-fc cases (right panel). We see the life-time of simulation 
becomes four-times longer than the standard ADM by tuning the parameter k. 



As one can check by (j46| and (|47|) adding Rij term 
keeps the constraint propagation in a first-order 
form. Frittelli (see also 76]) pointed out that 
the hyperbolicity of constraint propagation equa- 
tions is better in the standard ADM system. This 
stability feature is also confirmed numerically, and 
we set our CAF conjecture so as to satisfy this dif- 
ference. 

2. Detweiler type 

Detweiler [30| found that with a particular com- 
bination, the evolution of the energy norm of the 
constraints, H^+M^, can be negative definite when 
we apply the maximal slicing condition, K = 0. His 
adjustment can be written in our notation in 
and as 

= -KLa^-f,j, (49) 
= KLa\K,j - {l/3)Kj,,), (50) 
= «La2[3(a(,a)<5^)-(aza)7,,7''], (51) 

everything else is zero, where is a multiplier. De- 
tweiler's adjustment, (|15|) - ([5^ . does not put con- 
straint propagation equation to a first order form, 
so we cannot discuss hyperbolicity or the character- 
istic speed of the constraints. From perturbation 
analysis on Minkowskii and Schwarzschild space- 
time, we confirmed that Detweiler's system pro- 
vides better accuracy than the standard ADM, but 
only for small positive k^. 

We made various predictions how additional adjusted 
terms will change the constraint propagation [U Fm | . In 



that process, we applied CAF analysis for Schwarzschild 
spacetime, and searched that when and where the nega- 
tive real or non-zero imaginary eigenvalues of the homog- 
enized constraint propagation matrix appear, and how 
they depend on the choice of coordinate system and ad- 
justments. We found that there is a constraint-violating 
mode near the horizon for the standard ADM formu- 
lation, and that constraint-violating mode can be sup- 
pressed by adjusting equations and by choosing an ap- 
propriate gauge conditions. 



2. Numerical demonstrations and remarks 

Systematic numerical comparisons are in progress, and 
we show two sample plots here. Figure [7] is the case 
of Teukolsky wave [t^ propagation under 3-dimensional 
periodic boundary condition. Both the standard ADM 
system and Detweiler system [one of the adjusted ADM 
system with adjustments P5|) -([5 ^ ] are compared with 
the same numerical parameters. Plots are the L2 norm 
of the Hamiltonian constraint H.^^^^ , i.e. the violation of 
constraints, and we see the life-time of the standard ADM 
evolution ends up at f = 200. However, if we chose a 
particular value of [multiplier in we observe 

that violation of constraints is reduced than the standard 
ADM case, and simulation can continue longer than that 
(left panel). If we further tuned k^, say turn-off to = 
when the total L2 norm of Ji^^M jg gniall, then we can 
see that the constraint violation is somewhat maintained 
at a small level and more long-term stable simulation is 
available (right panel). 

During the comparisons of adjustments, we found that 
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FIG. 8: The one-dimensional gauge-wave test with the BSSN system (left) and the adjusted BSSN system (right) in the A- 
equation (|53p . The L2 norm of H, rescaled by the resolution parameter are plotted with a function of the crossing-time. 

The wave amplitude is set to 0.01, and we chose the adjustment parameter ka = 0.005. The BSSN system loses convergence at 
the early time, near the 20 crossing-time and it will produce the blow-ups of the calculation in the end, while in the adjusted 
version we see the higher resolution runs show convergence longer, i.e., the 300 crossing-time in Ti. and all runs can stably evolve 
up to the 1000 crossing-time. 



it is necessary to create time asymmetric structure of evo- 
lution equations in order to force the evolution on to the 
constraint surface. There are infinite ways of adjusting 
equations, but we found that if we follow the next guide- 
line, then such an adjustment will give us time asymmet- 
ric evolution. 

□ Trick to obtain asymptotically constrained system: 
= Break the time reversal symmetry (TRS) of the evolu- 
tion equation. 

1. Evaluate the parity of the evolution equation. 

By reversing the time {dt —dt), there are vari- 
ables which change their signatures (parity (— )) 
[e.g. Kij,daij,Mt,- ■ ■], while not (parity (+)) 
[e.g. g^j,dtKij,'H, ■ ■ ■]. 

2. Add adjustments which have different parity of that 
equation. 

For example, for the parity (— ) equation dtjij, add 
a parity (-I-) adjustment kTY. □ 

One of our criteria, the negative real CAFs, requires 
breaking the time-symmetric features of the original evo- 
lution equations. Such CAFs are obtained by adjusting 
the terms which break the TRS of the evolution equa- 
tions, and this is available even at the standard ADM 
system. 



C. Applications 2: Adjusted BSSN formulations 

1. Constraint propagation analysis of the BSSN equations 

In order to understand the stability property of the 
BSSN system, we studied the structure of the evolution 



equations, ([T7)l -(|2i p . in detail, especially how the mod- 
ifications using the constraints, (|22 p - (|^5|) . affect to the 
stability [t^I . We investigated the signature of the eigen- 
values of the constraint propagation equations, and ex- 
plained that the standard BSSN dynamical equations are 
balanced from the viewpoint of constrained propagations, 
including a clarification of the effect of the replacement 
using the momentum constraint equation, which was re- 
ported by Alcubierre et al. 

Moreover, we predicted that several combinations of 
modifications have a constraint-damping nature, and 
named them adjusted BSSN systems. Several adjusted 
BSSN systems are proposed in Table II of [77| . 

Yo et al. [ll] immediately applied one of our proposals 
to their simulations of stationary rotating black hole, and 
reported that one adjustment was contributed to main- 
tain their evolution of Kerr black hole [J/M up to 0.9M) 
for long time {t ~ 6000M). Their results also indicates 
that the evolved solution is closed to the exact one, that 
is, the constrained surface. 

Now, let us make clear some current technical tips 
listed in i jllB 21 using constraint propagation analysis. 

tip-1 Trace-out Aij technique can be explained that the 
violation of y^-constraint (PS)) affects to all other 
constraint violations. (See the full set of constraint 
propagation equations in Appendix of [z3-) 



tip-2 Replacement of T' enables to maintain tj-constraint 
Pi)l that delays the violation of U^SSN ^^^^ 
j^BSSN ^ (Again, the statement comes from the 
full set of constraint propagation equations. ) 
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2. Numerical Demonstrations 

We recently presented our numerical comparisons of 3 
kinds of adjusted BSSN formulation [4^. We performed 
three testbeds: gauge-wave, linear wave, and Gowdy- 
wave tests, proposed by the Mexico workshop on the 
formulation problem of the Einstein equations. We ob- 
served that the signature of the proposed Lagrange mul- 
tipliers are always right and the adjustments improve the 
convergence and stability of the simulations. When the 
original BSSN system already shows satisfactory good 
evolutions (e.g., linear wave test), the adjusted versions 
also coincide with those evolutions; while in some cases 
(e.g., gauge- wave or Gowdy-wave tests) the simulations 
using the adjusted systems last 10 times as long as those 
using the original BSSN equations. 

Figure [5] show the comparison between the (plain) 
BSSN system and the adjusted BSSN system in A- 
equation using the momentum constraint: 

dtA,j = dfA.j + KAaD(,Mj) , (53) 

where is predicted (from the eigenvalue analysis) 
to be positive in order to damp the constraint viola- 
tions. The testbed is one-dimensional gauge-wave; the 
trivial Minkowski space-time but sliced with the time- 
dependent 3-metric. The poor performance of the plain 
BSSN system of this test has been already reported [43 |. 
and one remedy is to apply the 4th-order finite differenc- 
ing scheme ^80] . The plots show that our adjusted system 
also improve the life-time of the plain BSSN simulation 
at least 10 times longer with better convergence. 



IV. OUTLOOK 

A. What we have achieved 

We reviewed recent efforts to the formulation problem 
of numerical relativity; the problem to find out a robust 
system against constraint violations. We categorized the 
approaches into 

(0) The standard ADM formulation C ^lTl^ . 

(1) The BSSN formulation (SHJ, 

(2) Hyperbolic formulations ( mi Cp . and 

(3) Asymptotically constrained formulations f ^II Dl) . 

Most of the numerical relativity groups now use the 
BSSN set of equations which are obtained empirically. A 
dramatic announcement of the success of binary black- 
hole simulations rushes the community to follow that 
recipe. Actually we are not yet completely understand- 
ing why the current set of BSSN equations together with 
particular combinations of gauge condition works well. 
Several explanations are applied based on the hyperbolic 



formulation scheme, but as we viewed there are not yet 
satisfactory. 

Our approach, on the other hand, tries to construct an 
evolution system that has its constraint surface as an at- 
tractor. Our unified view is to understand the evolution 
system by evaluating its constraint propagation. Espe- 
cially we proposed to analyze the constraint amplification 
factors which are the eigenvalues of the homogenized con- 
straint propagation equations. We analyzed the system 
based on our conjecture whether the constraint amplifi- 
cation factors suggest the constraint to decay/propagate 
or not. We concluded that 

• The constraint propagation features become differ- 
ent by simply adding constraint terms to the origi- 
nal evolution equations (we call this the adjustment 
of the evolution equations). 

• There is a constraint-violating mode in the stan- 
dard ADM evolution system when we apply it to 
a single non-rotating black hole space-time, and its 
growth rate is larger near the black-hole horizon. 

• Such a constraint-violating mode can be killed if 
we adjust the evolution equations with a particular 
modification using constraint terms. An effective 
guideline is to adjust terms as they break the time- 
reversal symmetry of the equations. 

• Our expectations are borne out in simple numeri- 
cal experiments using the Maxwell, Ashtekar, and 
ADM systems. However the modifications are not 
yet perfect to prevent non-linear growth of the con- 
straint violation. 

• We understand why the BSSN formulation works 
better than the ADM one in the limited case (per- 
turbation analysis in the fiat background), and 
further we proposed modified evolution equations 
along the lines of our previous procedure. Some 
of these proposed adjusted systems are numerically 
confirmed to work better than the standard BSSN 
system. 

The common key to the problem is how to adjust the 
evolution equations with constraints. Any adjusted sys- 
tems are mathematically equivalent if the constraints are 
completely satisfied, but this is not the case for numeri- 
cal simulations. Replacing terms with constraints is one 
of the normal steps when people re-formulate equations 
in a hyperbolic form. 

In summary, let me answer the following three ques- 
tions: 

• What is the guiding principle for selecting evolution 
equations for simulations in GR? 

-The key is to analyze the constraint propagation 
equation of the system. 

• Why many groups use the BSSN equations? 
-Because people just rush, not to be late to others. 
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• Are there an alternative formulation better than 
the BSSN? 

-Yes, there are. But we do not know which is the 
best one yet. 

B. Future directions 

If we say the final goal of this project is to find a robust 
evolution system against violation of constraints, then 
the recipe should be a combination of (a) formulations of 
the evolution equations, (b) choice of gauge conditions, 
(c) treatment of boundary conditions, and (d) numerical 
integration methods. We are now in the stages of solving 
these mixed puzzles. 

Recent attentions to higher dimensional space-time 
studies are waiting for numerical researches, but it is 
known that the formulation problem also exists in higher 
dimensional cases 66]. 

We have written this review from the viewpoint that 
the general relativity is a constrained dynamical system. 
This is not a proper problem in general relativity, but 
there is also in many physical systems such as electro- 



dynamics, magnetohydrodynamics, molecular dynamics, 
mechanical dynamics. Therefore, sharing and interact- 
ing the thoughts between different fields will definitely 
accelerate the progress. 

The ideal almighty algorithm to solve all the prob- 
lems may not exist, but the author believe that our final 
numerical recipe is somewhat an automatic system, and 
hope that numerical relativity turns to be an easy toolkit 
for everyone in near future. 
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